/*
 * FAtoms.cpp
 *
 *  Created on: May 13, 2012
 *      Author: marchi
 */

#include "FAtoms.h"
class brot{
	matrix co;
public:

	const matrix & operator()(const double a,const double b,const double c,
			const double alfa,const double beta,const double gamma){
	double degrad=M_PI/180.0;

	double ax=a;
    double alf=cos(degrad*alfa);
    double bet=cos(degrad*beta);
    double qt=sin(degrad*gamma);
    double gam=cos(degrad*gamma);
    double bx=b*gam;
    double by=b*qt;
    double cx=c*bet;
    double cy=c*(alf-bet*gam)/qt;
    double cz=sqrt(c*c-cx*cx-cy*cy);
    co[YY][XX]=0.0;
    co[ZZ][XX]=0.0;
    co[ZZ][YY]=0.0;
    co[XX][XX]=ax;
    co[XX][YY]=bx;
    co[XX][ZZ]=cx;
    co[YY][YY]=by;
    co[YY][ZZ]=cy;
    co[ZZ][ZZ]=cz;
    return co;
	}

};
struct duo{
	int n;
	double a;
};
struct op_duosort{
	bool operator()(const duo & x, const duo & y){return x.a < y.a;}
};
void FAtoms::ReadaStep(std::ifstream & fin){
}

void FAtoms::ReadaStep(FstreamC * fin){
	t_fileio *xd = fin->getfin();
	rvec   *x;
	matrix box;
	int    natos,step;
	real   prec,time;
	gmx_bool   bOK;
	int my_nframe=fin->gFrame();
	int nnframe=(my_nframe-fin->goffStep())/fin->gStep();
	try{
		if(firsttime){
			firsttime=false;
			if(my_nframe != fin->goffStep()){
				if((my_nframe-fin->goffStep())/fin->gStep() <= fin->gmFrameNumber()){
					xtc_seek_frame(fin->getfin(), my_nframe, nr);
				}
				else {
					cout << my_nframe << endl;
					throw string(" Beginning frame is out of range ");
				}
			}
		}
		natos=nr;
		x=new rvec[nr];
		if(!read_next_xtc(xd,natos,&step,&time,box,x,&prec,&bOK))
			throw string("Cannot read next frame!");
	}
	catch(const string & s){
		cout << s << endl;
		exit(-1);
	}
	Metric Met(box);
	setCoord(Met,x);
	doCOtoOC();
	fin->nextFrame();
	delete [] x;
}

void FAtoms::ReadaStep(FstreamF * fin){
	ReadaStepF(fin->getfin());
}

void FAtoms::ReadaStepF(std::ifstream & fin){
	double table[6];
	fin.seekg(FORTRANBYTES,ios::cur);
	fin.read(reinterpret_cast<char *> (&table),sizeof(table));
	fin.seekg(FORTRANBYTES,ios::cur);
	table[0]*=unit_nm;
	table[2]*=unit_nm;
	table[5]*=unit_nm;
	double alpha = 90.0 - asin(table[1]) * 90.0 / M_PI_2;
	double beta = 90.0 - asin(table[3]) * 90.0 / M_PI_2;
	double gamma = 90.0 - asin(table[4]) * 90.0 / M_PI_2;

	Metric Met(brot()(table[0],table[2],table[5],alpha,beta,gamma));
	vector<float> X[DIM];

	for(int i=0;i<DIM;i++){
		X[i].resize(nr,0.0);
		fin.seekg(FORTRANBYTES,ios::cur);
		fin.read(reinterpret_cast<char *> (&X[i][0]),sizeof(X[0][0])*nr);
		fin.seekg(FORTRANBYTES,ios::cur);
	}
	rvec * x=new rvec[nr];
	for(int n=0;n<nr;n++){
		for(int o=0;o<DIM;o++) x[n][o]=X[o][n]*unit_nm;
	}

	setCoord(Met,x);
	doCOtoOC();
}
void FAtoms::pdb(const vector<string> & data){

	vector<Dvect> xc;
	Metric Met(brot()(1000.0,1000.0,1000.0,90.0,90.0,90.0));
	for(size_t i=0;i<data.size();i++){
		if(data[i].find("ATOM") == 0 || data[i].find("HETATM") == 0 ) {
			Dvect tmp;
			std::stringstream ss(data[i].substr(30,24));
			ss>> tmp;
			xc.push_back(tmp);
		}
	}
	try{
		if(nr != xc.size()) throw string("Something wrong: No. atoms do not match");
	} catch(const string & s){
		cout << s <<endl;
		exit(1);
	}

	rvec * x=new rvec[nr];
	for(int n=0;n<nr;n++){
		for(int o=0;o<DIM;o++) x[n][o]=xc[n][o]*unit_nm;
	}
	setCoord(Met,x);
	doCOtoOC();
}

void FAtoms::moveOffset(FstreamC * fin){
	fin->nextFrame();
}
void FAtoms::moveOffset(FstreamF * fin){
	moveOffset(fin->getfin());
}
void FAtoms::moveOffset(std::ifstream & fin){
	const int DOUBLE=8;
	const int FLOAT=4;
	ios::streamoff OFFSET=FORTRANBYTES*2+6*DOUBLE+(FORTRANBYTES*2+nr*FLOAT)*DIM;
	fin.seekg(OFFSET,ios::cur);
}
void FAtoms::dodist(){
	int n=630-1;
	Dvect x0=x[n];
	vector<duo> d;

	for(size_t o=0;o< nr;o++){
		Dvect x1=x[o];
		duo tmp;
		tmp.n=o;
		tmp.a=x0.Dist(x1);
		d.push_back(tmp);
	}
	std::sort(d.begin(),d.end(),op_duosort());
	for(int n=0; n< 100; n++)
		cout << d[n].n+1 << " " << d[n].a*10.0 << endl;
}
FAtoms::~FAtoms() {
	// TODO Auto-generated destructor stub
}

